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ABSTRACT 

Late in their evolution, massive stars may undergo periods of violent instability and 
mass loss, but the mechanism responsible for these episodes has not been identified. We 
study one potential contributor: the development of local radiation-driven instabilities 
in the outer layers of main sequence (MS) and post-MS massive stars. We construct a 
sequence of massive stellar evolution models and investigate where they are subject to 
local radiative instabilities, both in the presence of magnetic fields and without them, 
and at a range of metallicities. We find that these types of instabilities do not occur in 
solar- metallicity MS stars up to 100 Mq, but they set in immediately post-MS for stars 
heavier than ~ 25 Mq. Once an instability appears, it involves a significant amount 
of mass in the star's upper layers (up to ~ Iper cent of the initial stellar mass), 
suggesting that radiation-driven instabilities are a potentially viable mechanism for 
dynamic mass loss. We find that the presence of magnetic fields at strengths low 
enough not to disturb the hydrostatic balance of the star does not alter these results. 
Stars with sub-Solar metallicity also show instability, but their instabilities involve 
less mass and appear later in the star's evolution. 

Key words: instabilities - stars: evolution - stars: interiors - stars: massive - stars: 
mass-loss - stars: magnetic field - supernovae: general 



INTRODUCTION 



The evolution of massive stars is driven lara 



jely by mass loss (e.g. IChiosi fc Maedeijll98d ). This loss takes the form of both 



continuously driven winds l|Puls et al.l 120081 . and references th erein) and episodes of violent mass ejection, althoug h which 
dominates the overall mass loss budget is still disputed (e.g. ISmith fc Owockil 120061 versus IVink fc Grafeneiil2012l ). While 
the continuo us winds can be understood as the produ ct of resonant scattering of photons off metal atoms in the stellar 
atmosphere I Castor et al.lll975l : IVink et al.ll200Gl . l200lh . the origin of the eruptive episodes is much less certain. Supernova 



observations, in particular of Type Iln, show evid ence for m assive st ars undergoing siza ble erup tive outbursts (> 10 ^M( 



within months to years prio r to their explosio n (e.g lDopita et al. 1984 : iChugai fc Danzigenil994 : Smith et al.ll2007l : lFolev et al. 
20071 : IPastorello et al.ll2012l : lOfek et ahlboish . 



T he prototypical example of violent mass loss is the Great Eruption of?; Car i nae in the 1840s l|Davidson fc Humphreys 
2997), which ejected several tens of Mq over a 20 year period ( Gomez et al.ll2006l . l2010|). The rate of ma ss loss rules out the 
line scattering that is responsible for the continuous winds as a possible mechanism (jOwocki et al.ll2004l 'l. Instead, proposed 
models fall into two broad classes. In one, the mass loss is powered by the internal radiation of the star interacting with 
a source of continuum rather than line opacity, and is thought to be triggered by the development of clumps in the stars ' 
upper layers that cau ses a rapid change in the effective continuum opacity (e.g. IShayiy||2000l . l2001bl : ISmith fc Owockil l2006l : 
van Marie et al.ll2008l l. The other class of models invokes tidal perturbation by a companion as the driving mechanism. The 
perturbation, applied to a star that is already on the edge of stability , induces an episode of mass transfer that gives rise to 
an outburst and mass loss (e.g. ISokeiil2004l . l2005l : iKashi fc Sokeill20ld ). 
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Either of the proposed mechanisms calls for the star to be on the verge of instability, so that some perturbation, either 
internal or external, can trigger an outburst. Since the amount of mass involved is large, this instability cannot be confined to 
the stellar atmosphere but must involve deeper layers o f the star. There ar e several candidate global instabil ities that affec t 
stars near the Eddington lim it, including strange modes |Glatzel et al.lll999l ). global radiatively-driven modes ( Shavivll2001al ). 
and opacity-driven g-modes (JTownsend fc MacDonaldl 120061 ). 

However, thus far there has been no significant search for local instabilities that might play a role in destabilizing the 
outer layers of massive stars, making them ripe for eruptive ejection. This is somewhat surprising, because the criterion for 
local radiation-hydrodyn amic and radiation - magn etohydro dynamic instabilit ies in optically thick, gravity-confined media was 
derived very generally bv lBlaes &: Socrated (|2003). Indeed, Blaes fc Socrated conjecture that the outer layers of massive stars 
might be subject to their computed instability I Fernandez fc Socrateal2012l ). 

In this paper we investigate that conjecture, and determine where, in terms of both stellar mass and evolutionary state, 
local radiation-driven instabilities in stars will occur. We construct a grid of evolutionary tracks for the interiors of stars 
at a variety of masses using the MESA (Modul es for Experiments in St ellar Astrophysics) package IjPaxton et aLll2011f ). and 
determine where in that grid stars satisfy the lBlaes fc Socrated (2003|) instability condition. We show that local instability 
does not occur in main sequence stars, but that it appears almost immediately after the end of hydrogen burning in stars with 
initial masses above ^ 25 Mq. When present, the instability occurs at high optical depth within the star, and the unstable 
region typically contains up to 1 Mq of material, suggesting that this instability is a candidate contributor to eruptive mass 
loss. 

The remainder of this paper is structured as follows. In Section [2] we present the details of our models and stability 
calculations. We report the results of this exercise in Section [S] and discuss their implications in Section [l] 



2 METHODOLOGY 

In this section, we first describe the instability criteria we will use, and then discuss how we generate the stellar structure 
models to which we will apply them. 



2.1 Radiation-driven hydro- and magneto-hydrodynamical instabilities 

Starting from the equations of radiation hydrodynamics in the difi'usion limit, iBlaes fc Socrated ( 20031 ) obtain dispersion 
relations for local radiation pressure- driven instabilities in the local Wentzel-Kramers-Brillouin (WKB) limit. Their analysis 
relies on a few assumptions: in the unperturbed state (but not necessarily in the perturbed state) matter and radiation are 
in LTE, the matter is of constant composition and ionization state, gravitational potential perturbations are negligible (the 
Cowling approximation), as is photon viscosity. They consider both purely hydrodynamic (HD) and magneto-hydrodynamic 
(MHD) cases; for the latter, the magnetic field in the unperturbed state is required to be constant. 

For the HD case, the dispersion relation for the acoustic waves (w oc fc), to first order, is given by equation (62) of 
Blaes fc Socrated ( 20031 ) and reproduced here: 
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where k is the wave vector, Ci is the isothermal sound speed in the gas, kf is the flux mean opacity, c is the speed of light, 
Pgas is the gas pressure, E and F are, respectively, the radiation energy density and fiux, and 0p is the logarithmic derivative 
of K,F with respect to the density, i.e. 
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From this dispersion relation, the condition for instability is 
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For t; > 0, t; is the growth rate of the instability. Since all physical quantities involved are positive, instability occurs only when 
the second term inside square brackets is larger than the first term and opposite in sign, that is, for the upward propagating 
acoustic waves corresponding to the minus in the plus-minus sign. The first term represents damping by radiative difi'usion, 
while the second represents forcing by the radiation fiux. Note that the forcing term is proportional to the logarithmic 

derivative of the opacity Op. This will become significant in our analysis below. 

The dispersion relation for magneto-acoustic waves, to first order, is given by equation (92) of iBlaes fc Socrated (2003|) 
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and reproduced here: 
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where va is the magnitude of the vector Alfven speed (va = B/[47rp] ' ), and Wph is the phase speed of the wave, defined as 
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This dispersion relation gives the condition for instabihties to occur, i.e. 
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again, when instabihty occurs, r\ gives the growth rate. We note that in this case there are two terms involved in the condition: 
the first one is a magnetically-modified hydrodynamic term, analogous to the one presented in equation ([3]) and in which the 
opacity driving must overcome the radiative diffusion to contribute to an instability. However, the inclusion of a magnetic 
field in the analysis introduces the factor 
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which is always positive. The second term in equation ((S)) corresponds to a purely magnetohydrodynamic contribution, and 
is positive both for "fast" and "slow" magnetosonic waves. As such, this term always favors instability, but will be important 
only when it dominates the hydrodynamic term. 



2.2 Stellar models in MESA star 

We will use the instability conditions given by equations ([Sjl and (|6} on stellar profiles of stars computed using MESA star to 
determine whether instabilities might be able to develop in stars at a range of masses and ag es. MESA star is a o ne-dimensional 
stellar evolution code included in the MESA package for computational stellar astrophysics (|Paxton et al.ll201ir ). MESA star is 
a Henyey sty l e cod e that simultaneously solves the fully coupled structure and composition equations; we refer readers to 
Paxton et al.l I 201l| ) for details on its capabilities and limitations. 

We use MESA star to simulate 21 stars with different masses, all with Solar metalicity, using the default code parameters. 
Our model stars have masses from SM© to 100 M©, in intervals of 5 Mq, plus a model for a IM© star. This range in masses 
allows us to study differences in behavior between low mass and high mass stars. In order to understand how metallicity 
affects the results, we also perform a similar set of runs at metallicities of 1, 3, 10, and 20 per cent of Solar metallicity. At 
each output time for each mass, we evaluate the instability conditions given by equations ([3]) and ((6]) at every position within 
the structure of the calculated stellar profiles. We take the partial derivatives of the opacity Qp directly from the Mesa star 
output. 

There are two complications to our analysis. First, MESA star does not include magnetic field s. To evaluate the MHD 
instability condition, we therefore adopt constant magnetic fields at a range of strengths. Second. iBlaes fc Socrated ( 20031 ) 
only derive their instability criteria for conditions of constant mean molecular mass, and this assumption breaks down within 
stars' hydrogen and helium ionization zones. For this reason, we exclude from our analysis all those locations where the 
hydrogen ionization fraction is 1 — 99 per cent, or where the helium ionization fraction is 10 — 90 per cent. This ensures that 
we are only considering regions where the mean molecular mass is constant to 1 per cent or better. Note that this exclusion 
does not mean that there cannot be instabilities present in the ionization zones, simply that the instability criteria to which 
we have access are not applicable in those regions. 



3 RESULTS 

3.1 Hydrodynamic Instabilities in Single Stars 

We first consider the purely HD case at Solar metallicity; we defer discussion of magnetic and metallicity effects to Section [3]3l 
Figure [l] shows sample results from one of our models: we plot the physical properties of a 70 M0 star as a function of the 
mass coordinate at two different moments of its life. The quantities shown are those required to evaluate the HD instability 
condition (equation[3|. The upper panel shows the state of the star while it is on the main sequence; as the plot shows, at this 
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M = 70Mg 

age=2.931Myr 

MS lifetime=3.156Myr 




64 66 



Figure 1. Properties of a star with an initial mass of 7OM0as a function of the mass coordinate in two stages of its evolution. All 
dimensional quantities are measured in cgs units. The upper panel shows the state of the star at an age of 2.931 Myr, while the star is 
still on the main sequence, while the lower panel shows its state at 3.477 Myr, after the end of hydrogen burning at 3.156 Myr. Note 
that the mass of the star is < 70 Mq due to stellar winds. The quantities plotted are those required to evaluate the HD growth rate rj 
(equation [3J , including rj itself. Dashed black lines indicate regions where ?; < and the star is stable, while solid black lines indicate 
regions where jy > and the star is unstable. The gray shaded bars mark the hydrogen and helium ionization zones that we exclude from 
our analysis; see Section l2.2l The upper panel shows a star without instabilities, while the lower one presents a stage in which the star 
has already developed an opacity-driven instability. 



time the system does not yet show any relevant instability. As discussed above, the condition for an instability to be present 
is that the radiative flux term must overcome diffusive damping; the former is proportional to the logarithmic derivative of 
the opacity with density Qp, and during the main sequence stage Op ~ 0.01 throughout most of the star. The underlying 
physical reason for this small value is that, on the main sequence, the stellar material is very hot even fairly far out in the star. 
As a result, the opacity is dominated by electron scattering which is density-independent, rather than free- free or bound- free 
processes which are not. As a result, radiative forcing is weak and the star is stable. 

The lower panel, on the other hand, shows the star after it leaves the main sequence. At this point a n acoustic instability 
has developed, even excluding the ionization zones where we cannot apply the iBlaes fc Socrated ( 20031 ) instability analysis. 
Instability occurs in two zones near the surface of the star, on either side of the hydrogen and helium ionization zones. There 
is a notable difference between the unstable zones: in the innermost, the radiation pressure E and the gas pressure Pgas are 
at the same order of magnitude, while in the region that reaches the stellar surface Pgas is much greater than E. The total 
mass contained in the two unstable zones together is ~ 1 M0 at this evolutionary stage (Figure [!}. The underlying reason for 
the instability is that, as the star leaves the main sequence, it expands and the gas temperature in the outer layers drops 
drastically; this causes free-free and bound-free opacity to become competitive with electron scattering, which in turn raises 
Op and destabilizes the gas. 
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Figure 2. Growth rate (tj) for HD instabilities in the outer regions of a 70 Mq star, as a function of its age. Gray regions indicate where 
r; < 0, so the star does not satisfy the HD instabiUty condition. Color regions indicate values of 77 > 0. Note that the vertical axis has 
three different linear scales for the age (separated by white dashed lines) in order to show the behavior of different evolutionary stages of 
the star: its MS (up to 3.156Myr), a post-MS stage in which the star gradually loses mass and builds a layer of unstable material (from 
3.156Myr to 3.475Myr) and a later stage in which the star becomes more unstable towards its surface (from 3.475Myr to 3.478Myr). 



3.2 Hydrodynamic Instabilities as a Function of Stellar Mass and Age 

By repeating the analysis above for a given star at different stages of its evolution, we can find where and when the instability 
appears for a given initial mass. Figure [2] shows the growth rate of the HD instability computed over the outer 15 per cent 
of the mass of our initial 70 Mq star as a function of its age. Negative values for rj, indicating stability, are shown in grey. 
As the star loses mass, a white region appears on the right side of the figure, representing this decrement. The end of the 
main sequence for the star is shown with the lower white dashed line and a label. To accommodate in a single graph the 
varying timescales present during the stellar evolution, the age axis is broken in three different linear scales, separated by 
the two white dashed lines in the image. This figure shows that the HD instability condition (equation |3} is fulfilled in the 
outermost region of the star almost immediately after it finishes core hydrogen burning. The amount of mass contained in the 
unstable region tends to grow as the star evolves. The inverse of the growth rate is a characteristic time of the instability; for 
most of the unstable zone, the growth time is smaller than ~ 10* s (i.e. rj > 10~* s~^), which amounts to a few years and is 
significantly less than the age of the star. We also find that the wavelength of the most unstable mode is much smaller than 
either the stellar radius or the radial extent of the unstable region, indicating that instabilities will not be quenched by the 
finite thickness of the unstable layer. 

By extending our analysis to stars of different masses, we can compute the mass Mu in the unstable region for each star 
as a function of initial mass Afo and stellar age. Note that our calculation of AIu excludes the ionization zones, so it provides 
a lower limit for the size of the unstable region. Figure [3] shows our results. For both panels, the horizontal axis shows initial 
mass, while the vertical axis represents the age of each star. In the upper panel we measure age in units of each model's MS 
lifetime, defined as the time at which hydrogen burning in the core ceases; the white dashed horizontal line indicates this 
moment (age = agejvis)- The lower panel shows the same result, but now with the vertical axis indicating absolute rather 
than normalized age. The upper limit of the vertical axis in this plot is ~ 14 Gyr, comparable to the age of the Universe; we 
therefore see that a 1 Mq star will not become unstable over the entire age of the Universe. From these maps, it is clear that 
although most of the stars reach a stage where they have appropriate conditions for acoustic radiative-driven instabilities, it 
is only massive stars (above ~ 25 M0 ) that actually have a significant fraction of their mass contained in the unstable region, 
and where the instability appears rapidly after the end of the main sequence rather than extremely late in their lives. 

It is also of interest to investigate how deep within the star the unstable regions lie. In Figure |4] we plot the Rosseland 
mean optical depth r from the stellar surface to the inner edge of the outermost unstable region. The coordinate axes in the 
figure are the same as those in the upper panel of Figure |3l age in units of MS age versus initial stellar mass. This graph 
demonstrates that the unstable regions of interest have a high optical depth, which suggests that this mechanism can produce 
a layer of unstable mass lying deeper than the stellar photosphere at r = 2/3. It is interesting to note that even the models 
in which the fraction of mass in the outer unstable region is minuscule, the optical depth is quite high. When instabilities are 
present, they typically go down to r ~ 10. 
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Figure 3. Mass contained in the unstable regions (Afu) of stars as a function of their initial mass (horizontal axis) and age (vertical 
axis). Upper panel: the age of each star is shown in units of its MS lifetime. Lower panel: ages are shown in years, with the highest value 
plotted corresponding to the current age of the Universe. 
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Figure 4. Colors show the Rosseland mean optical depth (t) from the stellar surface to the inner edge of the outermost unstable region. 
Axes are identical to those in the upper panel of Figure |3] Unstable regions extend significantly deeper than the stellar atmosphere 
(defined at t = 2/3). 

3.3 Effects of Metallicity and Magnetic Fields 

All the stellar models we have considered so far have Solar metallicity: Zq. Here we extend our formalism to include sub-Solar 
metallicity models with 0.2, 0.1, 0.03, and 0.01 Zq. We find that the qualitative result that the mass in the unstable region 
increases with initial stellar mass, while the time post-MS at which instability appears decreases, continues to hold for varying 
metallicity. However, the amount of mass in the unstable region (again excluding the ionization zones) decreases for all models 
as the metallicity decreases. We illustrate this effect in the upper panel in Figure O which shows a map similar to that plotted 
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Figure 5. Color maps showing tlie amount of mass contained in the outermost unstable regions of each model at different stages of 
their evolution. The upper panel shows models for which the metallicity has been decreased to 0.03 Zq: the amount of unstable mass 
involved has decreased in all stars, but the general tendency seen before remains. The models depicted in the lower panel are identical 
to the ones in the upper panel in Fig.|3l but a constant-magnitude magnetic field of 100 G has been included. The effect in the unstable 
masses is indistinguishable. 



in the upper panel of Figure [3l but for 0.03 Z©. The general weakening of the instability at lower metallicity can again be 
understood by considering the density-dependence of the dominant opacities, and thus the value of 0p. As the metallicity 
decreases, the supply of heavy ions in the outer layers of the star decreases, and thus the importance of bound-free opacity is 
reduced compared to the electron scattering opacity. Since the former is density-dependent and the latter is not, the overall 
opacity becomes less density-dependent (0p is reduced) and the radiative force is diminished. This weakens the instability. 

We also investigated the MHD instability condition (equation [6l ) using constant-magnitu de magnetic fields ranging from 
1 G to 10^ G. In order to maintain the static equilibrium presumed bv lBlaes fc Socrated (2003), and to remain within the range 
of validity of the non-magnetic MESA star models, we only consider magnetic fields small enough that the magnetic pressure 
is significantly smaller than the gas pressure. A lO'' G magnetic field amounts to a pressure of « 4xl0^''Ba, significantly 
lower than the gas and radiation pressures present during the MS of a 70 Mq star, as seen in the upper panel of Figure [T] on 
the other hand, the strongest magnetic field that we could safely propose for the same star when the instabilities are present 
is of the order of 100 G, corresponding to a magnetic pressure of « 400 Ba (see lower panel of Figure [!}. This restriction is 
reasonable, since there is no evidence that magnetic fields in stellar interiors are ever strong enough for magnetic pressure 
to significantly alter their hydrostatic balance, and it seems unlikely that one could construct a hydrostatic spherical stellar 
model if they were. With this restriction in mind, we find that magnetic fields in the allowed range of strengths for already 
developed instabilities have negligible effects on the size of the instability region. The hydrodynamical contribution to the 
instability always dominates. The lower panel in Figure [5] presents the case where a magnetic field of magnitude B — 100 G 
was set up. This map is practically indistinguishable from its HD counterpart (upper panel in Fig. [3]) and illustrates the lack 
of a magnetic effect in this situation. Applying stronger fields (up to 10® G) to a 70 M0 star during its MS does not contribute 
to the formation of instabilities. 

We can understand these results by examining the form of the instability condition. Late in the star's evolution, when 
radiative forcing dominates the hydrodynamic term in equation ((6|, the ratio between the hydrodynamic and magnetic terms 



is of order (vph/vA) 6p. If the magnetic field is dynamically weak compared to gas pressure, then iiph ~ Ci and the term 
becomes of order {ci/vA)'^Op- Thus the magnetic contribution to the instability condition is significant only if 6p ^ (vA/ci)'^. 
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Figure [1] shows that Qp is typically of order a few tenths in the outer parts of post-MS stars. In this post-MS case, this 
means that there is no room for magnetic fields to have a significant effect. The requirement that magnetic pressure be small 
compared to gas pressure implies that (vAJCif" ^ 1, but magnetic fields only significantly alter the instability condition if 
{vA/ci)^ ^ Qp ~ 0.1. This analysis shows that it is clearly impossible for the magnetic field to have anything more than a 
marginal effect. In the MS case, where the HD term is dominated by radiative diffusion damping, the ratio of HD and magnetic 
terms is of order (cf/t;^)ci(i5 + Pgas)/-?, assuming again that the magnetic field is weak compared to the radiation and gas 
pressures. In this case we can make a similar argument. The radiative flux F is of order cE/t, where r is the characteristic 
optical depth, which means that the factor Ci{E + Pgas)/-F ~ T{Pgas/E){ci/c). The factor Ci/c ~ 10~^ is small, but deep in 
the stellar interior r ^ 10"^, and far out in the star Pgas/-B S> 10^, so the total factor is everywhere ^ 1; examination of 
Figure [T] indicates that it is of order 10* — 10®. Th is large value is a dire ct consequence of the fact that stellar interiors are in 



the dynamic rather than the static diffusion limit (jKrumholz et al.ll2007f ). At the same time, the restriction that the magnetic 



field not dominate the hydrostatic balance of the star implies that (cf/w^) > 1. Together these two constraints imply that 
the HD term must greatly exceed the magnetic one in main sequence stars. In practice and in theory, we find that including 
magnetic fields has little practical effect in the behavior of the instabilities already present in the purely HD case. 



4 CONCLUSIONS AND IMPLICATIONS 

In this paper we investigate under what conditions local radiation-driven instabilities will be present i n stellar interiors. T o 
this end we construct a series of stellar models as a function of mass and age using the MES A star code JPaxton et al.ll201ll l 



and then apply the hydrodynamic and magnetohydrodynamic instability criteria derived bv lBlaes &c Socrated (2003). We find 



that that the outer layers of massive stars can present the properties required for instability to develop. This does not occur 
in main sequence stars up to 100 M0, but it sets in immediately post-MS for stars more massive than ~ 25 M©. Once these 
instabilities develop, they involve a non-negligible amount of the stellar mass, reaching up to ~ 1 M0 in the more massive cases 
studied. They also occur relatively deep in the stellar interior, at optical depths of order r — 10. The underlying physical cause 
for the onset of the instability is a drop in gas temperature as stars expand in the post-MS phase, which raises the importance 
of bound-free and free-free opacity compared to electron scattering. Bound-free and free-free opacity are density-dependent, 
and this density-dependence is destabilizing. 

This process could be a potential mechanism for the enhanced dynamic mass loss that is observed in some massive stars 
during their lifetimes and prior to exploding. The presence of a layer of material that is unstable in this way could prepare the 
conditions of an episode of eruptive mass loss, perhaps triggered by some external event such as the gravitational interaction 
of the star with a potential binary companion. At a minimum, the presence of the instability will modify the structure of 
the star in the regions where it occurs. However, determining the full results of the instability will require detailed numerical 
simulations. 
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